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ABSTRACT 

Several scenarios have been proposed to explain the presence of multiple stellar populations 
in globular clusters. Many of them invoke multiple generations of stars to explain the ob¬ 
served chemical abundance anomalies, but it has also been suggested that self-enrichment 
could occur via accretion of ejecta from massive stars onto the circumstellar disc of low-mass 
pre-main sequence stars. These scenarios imply different initial conditions for the kinematics 
of the various stellar populations. Given some net angular momentum initially, models for 
which a second generation forms from gas that collects in a cooling flow into the core of 
the cluster predict an initially larger rotational amplitude for the polluted stars compared to 
the pristine stars. This is opposite to what is expected from the accretion model, where the 
polluted stars are the ones crossing the core and are on preferentially radial (low-angular mo¬ 
mentum) orbits, such that their rotational amplitude is lower. Here we present the results of 
a suite of A^-body simulations with initial conditions chosen to capture the distinct kinematic 
properties of these pollution scenarios. We show that initial differences in the kinematics of 
polluted and pristine stars can survive to the present epoch in the outer parts of a large frac¬ 
tion of Galactic globular clusters. The differential rotation of pristine and polluted stars is 
identified as a unique kinematic signature that could allow us to distinguish between various 
scenarios, while other kinematic imprints are generally very similar from one scenario to the 
other. 

Key words: galaxies: star clusters: general - globular clusters: general - stars: kinematics 
and dynamics 


1 INTRODUCTION 

Once considered as quintessential simple stellar populations (i.e. no 
dispersion in age or metal content), globular clusters (GCs) now ap¬ 
pear to host multiple stellar populations. This has been revealed by 
ubiquitous star-to-star abundance variations of light elements such 
as Na, O, and A1 (e.g. Gratton, Carretta & Bragaglia 2012), indicat¬ 
ing that a significant fraction of cluster stars (~ 30 - 70%; e.g. Car¬ 
retta et al. 2009b) must have been polluted by material processed 
in high-temperature stellar interiors via proton capture reactions. 
It is also apparent from the multiple or extended main sequences, 
subgiant and red-giant branches seen in the colour-magnitude di¬ 
agrams of a growing number of Galactic GCs (e.g. Piotto et al. 
2012), which have been attributed to different levels of helium en¬ 
richment (e.g. D’Antona et al. 2002; Norris 2004; Pasquini et al. 
2011) and light element abundance variations (e.g. Marino et al. 
2008; Sbordone et al. 2011) among the cluster stars. 

Several models have been proposed in an attempt to explain 
the observed chemical abundance anomalies, with different types 
of stars suggested as the source of helium-enriched material. The 
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yields from these potential enrichment sources must also explain 
the observed chemical abundance patterns (e.g. the Na-O anticor¬ 
relation) without producing a spread in iron abundances (such a 
spread is absent in almost all but the most massive GCs; e.g. Car¬ 
retta et al. 2009a). Asymptotic Giant Branch (AGB) stars (e.g. Ven¬ 
tura et al. 2001; D’Ercole et al. 2008), fast rotating massive stars 
(“spin-stars”; e.g. Decressin et al. 2007), interacting massive bina¬ 
ries (de Mink et al. 2009) and super-massive stars with M ~ 10"* Mg 
(Denissenkov & Hartwick 2014) have all been considered as possi¬ 
ble polluters. 

Obtaining the right chemistry for the polluting material is one 
aspect of the problem, but finding a way to get this material into 
stars is another equally important and challenging issue. Many of 
the models put forward invoke the formation of a second generation 
of stars from the processed material released by a first generation. 
For this reason, the terms multiple generations and multiple pop¬ 
ulations are often interchangeably used in the literature to refer to 
the presence of polluted and pristine stars, but note that the notion 
of multiple star-formation events within a single cluster is an as¬ 
sumption of these models for which there is no direct observational 
evidence. 

In a scenario that has received significant attention, D’Ercole 
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et al. (2008) explored the possibility that a second generation is 
formed from the gas ejected by AGB stars. They used ID hydro- 
dynamical simulations to show that the low-velocity ejecta of AGB 
stars can be retained by the cluster and form a cooling flow that 
rapidly collects in the innermost regions. The new generation then 
forms a more concentrated stellar subsystem, in keeping with ob¬ 
servations of a more centrally concentrated polluted (i.e. Na-rich, 
O-poor) population in a number of Galactic GCs (e.g. Lardo et al. 
2011). This scenario was further studied by different authors (e.g. 
Conroy & Spergel 2011; Bekki 2011), many of which stressed the 
importance of accreting pristine interstellar gas onto the cluster for 
it to work. In particular, it was suggested that the total mass of ac¬ 
creted pristine gas must be comparable to that of AGB ejecta in or¬ 
der to reproduce the observed Na-O anticorrelation (e.g. D’Ercole 
et al. 2010; Ventura et al. 2013). How such pristine interstellar ma¬ 
terial can remain free of pollution by supernova ejecta from the first 
generation (as required by the lack of an iron abundance spread) 
however remains unclear. In addition to the above model based on 
AGB stars, the main scenario invoking spin-stars as the source of 
pollution also works under the assumption that the slow ejecta from 
these stars form the basis of a new stellar generation (e.g. Decressin 
et al. 2007). Charbonnel et al. (2014) even postulated that, within 
this framework, all low-mass stars present in globular clusters to¬ 
day could in fact be second-generation stars. 

The models that require multiple star-formation events come 
with many potential shortcomings (for an overview see Renzini 
2008; Bastian et al. 2013b; Cabrera-Ziri et al. 2015). For exam¬ 
ple, they generally suffer from an “internal mass budgef’ problem, 
as the amount of processed material returned by either AGB stars 
or spin-stars is insufficient to account for the large observed frac¬ 
tions of polluted stars in GCs. To circumvent this issue, these mod¬ 
els either invoke a non-standard initial mass function for the first 
or second generation (e.g. Decressin et al. 2007) and a 100% star- 
formation efficiency for the second-generation stars (e.g. D’Antona 
et al. 2013), or imply that GCs were initially ~ 10- 100 times more 
massive than observed at the present time and that a large frac¬ 
tion of first-generation (i.e. pristine) stars have since been lost (e.g. 
D’Ercole et al. 2008; Schaerer & Charbonnel 2011; Bekki 2011). 
The latter possibility appears to be in conflict with the “external 
mass budget” constraints imposed by observations of GCs and field 
stars in nearby dwarf galaxies (Fornax, WLM, IKN), which imply 
that GCs in these systems could not have been more than about 
five times more massive initially (Larsen, Strader & Brodie 2012; 
Larsen et al. 2014). These models are also at odds with observa¬ 
tions of young massive clusters forming at the present epoch, for 
which there is currently no conclusive evidence for age spreads, 
ongoing star formation, or reservoirs of gas and dust from which a 
new generation could form (e.g. Kudryavtseva et al. 2012; Bastian 
& Silva-Villa 2013; Bastian et al. 2013a; Cabrera-Ziri et al. 2014; 
Bastian & Strader 2014; Bastian, Hollyhead & Cabrera-Ziri 2014). 

To address these problems, Bastian et al. (2013b) proposed 
a new scenario in which processed material shed from interacting 
massive binaries during the first ~10 Myr of the cluster’s life is 
accreted onto the circumstellar discs of low-mass stars (still on the 
pre-main sequence), ultimately polluting those stars. It must be said 
that the suggestion that pollution could originate from the accre¬ 
tion of ejecta from higher mass stars was first made by D’Antona, 
Gratton & Chieffi (1983), who at the time considered that the ob¬ 
served abundance patterns were reflecting surface contamination. 
More recent evidence however indicates that abundance anomalies 
are present throughout the star, as evolved red giant branch stars 
(mixed through convection) show the same Na-0 anti-correlations 


as main-sequence stars (Cohen, Briley & Stetson 2002). The varia¬ 
tion proposed by Bastian et al. (2013b) addresses this by polluting 
stars when they are on the pre-main sequence and convective, and it 
also maximises the efficiency of accretion by invoking gas capture 
by discs. 

In the so-called “early disc accretion” scenario, the polluted 
stars originate from the same generation as the pristine stars, with 
no need for multiple star-formation events. Only the stars that pass 
through the cluster core to sweep up gas are polluted, as this is 
where high-mass stars and their ejecta are expected to be preferen¬ 
tially located, either because they form there (e.g. Bonnell, Bate & 
Zinnecker 1998) or because they segregate soon after their forma¬ 
tion (e.g. Portegies Zwart, McMillan & Gieles 2010). Depending 
slightly on the velocity structure and the gas distribution, this leaves 
about half of the low-mass stars, which did not cross the core, with 
normal abundances. The model also predicts an enriched popula¬ 
tion that is more centrally concentrated than the pristine population. 
Given that a large fraction of high-mass stars are in binary systems 
that will interact at some point in their evolution (Sana et al. 2012, 
2013), these interacting systems are promising sources of polluting 
material. It has been argued that their ejecta can display the required 
chemical properties and also satisfy the internal mass budget con¬ 
straints (de Mink et al. 2009) without requiring GCs to have been 
significantly more massive initially. 

Among the caveats of the early disc accretion model, an im¬ 
portant one is the uncertainty on the expected lifetime of circum¬ 
stellar discs around low-mass pre-main sequence stars in the dense 
environment of a globular cluster. These discs need to survive for 
5-10 Myr to sweep up enough material for this scenario to work, 
and Bondi-Hoyle accretion is expected to be inefficient due to 
the large velocity dispersions of massive clusters. D’Antona et al. 
(2014) also pointed out that the structure of the seed stars may 
not remain fully convective for the whole duration of the accre¬ 
tion phase. If this is true, the full mixing of the helium-enriched 
material required down to the stellar centre would not be achieved. 
Thermohaline mixing may offer a way around this, but the same 
authors suggested that this could destroy any lithium surviving in 
the envelope, which may be in contradiction with observations (e.g. 
Shen et al. 2010). 

Both the better-studied multiple generations model (with AGB 
stars as the source of pollution; e.g. D’Ercole et al. 2008) and the 
newer early disc accretion model (Bastian et al. 2013b) have their 
difficulties, but neither has been conclusively ruled out to date. In 
the present paper, we therefore compare two main scenarios for 
how enriched material makes its way into stars: one in which a sec¬ 
ond generation is formed from gas that collects at the centre of the 
cluster (hereafter referred to as the multiple generations scenario), 
and another one in which pollution is due to accretion onto low- 
mass stars from the same generation when they pass through the 
core of the cluster (hereafter referred to as the accretion scenario). 
Instead of focusing on the chemistry as is often the case in studies 
of multiple populations, we approach this problem from another 
angle and explore the fossil kinematic imprints that may follow 
from these two formation scenarios after a Hubble time of dynam¬ 
ical evolution. Our results are not tied to a specific source of pol¬ 
luting material. They apply to any “multiple generations” scenario 
that proceeds as suggested by D’Ercole et al. (2008), whether AGB 
stars are the polluters or not. Similarly, they would be relevant for 
any “accretion” scenario that proceeds as outlined by Bastian et al. 
(2013b), regardless of the source of pollution or whether accretion 
proceeds via discs or not. 

In §2, we start by reviewing previous work on the dynamics of 
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multiple populations, which was mainly set in the context of mod¬ 
els invoking multiple generations. §3 describes the Al-body simula¬ 
tions that we performed and the different initial conditions that we 
chose to capture the distinct kinematic properties implied by the 
two self-enrichment models considered. We then present in §4 the 
results of these simulations, with an emphasis on the identification 
of a unique kinematic signature that allows to distinguish between 
the two models. In §5, we discuss how this kinematic signature 
could be observed and we identify a large subsample of Galactic 
GCs in which it is expected to have survived to the present epoch. 
We present our conclusions in §6. 


2 PREVIOUS STUDIES OE THE DYNAMICS OE 
MULTIPLE POPULATIONS 

As mentioned previously, the internal mass budget problem of 
models invoking multiple generations can in principle be avoided 
by assuming that GCs were initially much more massive. A signifi¬ 
cant fraction of first-generation stars (~ 90% or more) must then be 
lost (i.e. unbound) by dynamical evolution while the more concen¬ 
trated second generation remains relatively unaffected. Decressin, 
Baumgardt & Kroupa (2008) showed with A-body simulations that 
two-body relaxation alone cannot cause such a strong preferential 
loss of first-generation stars, and suggested that a more efficient 
mechanism such as primordial gas expulsion is needed to unbind 
the first-generation stars located in the outer regions of the cluster 
on a short timescale. 

Decressin et al. (2010) found that, under favourable circum¬ 
stances, the expulsion of gas left over from star formation can lead 
to a cluster containing 60% of second-generation stars. The radial 
distribution of the two populations is still expected to be different 
after gas expulsion, so two-body relaxation could further increase 
the fraction of second-generation stars. Note however that the im¬ 
portance of early gas expulsion in young massive clusters has re¬ 
cently been questioned by hydrodynamical simulations (Kruijssen 
et al. 2012) and kinematic observations (e.g. Henault-Brunet et al. 
2012a; Cottaar et al. 2012) suggesting that newly formed clusters 
are relatively gas-poor at their centre and in virial equilibrium from 
a young age (a few Myr). 

D’Ercole et al. (2008) studied the effect of mass loss from 
type II supernova ejecta (without including primordial gas left over 
from star formation) and suggested that this alone could trigger 
rapid early expansion of the cluster and lead to a strong preferen¬ 
tial loss of first-generation stars. In any case, it seems that the initial 
conditions and mass lost in this early phase need to be fine-tuned 
to explain that the mass ratio of the observed populations is always 
of order unity. One may otherwise expect GCs experiencing weak 
tidal fields to contain a much larger fraction of pristine stars. This 
does not seem to be the case for the remote and massive cluster 
NGC 2419 (Ra = 94.7 kpc. My = -9.5). Multi-band photometric 
observations of this cluster indeed suggest that a significant fraction 
of its stars belong to a helium-enriched population (Beccari et al. 
2013). 

In the following sections, we will nevertheless ignore these 
caveats about early mass loss and follow the long-term dynamical 
evolution of multiple populations starting from a cluster in virial 
equilibrium with a similar number of polluted and pristine stars. 
This approach was also adopted by Vesperini et al. (2013), who 
studied the long-term evolution of the relative spatial distribution 
of first- and second-generation stars (with spherically symmetric 
distributions) in the context of the scenario presented by D’Ercole 


et al. (2008). Starting with a second generation that is more cen¬ 
trally concentrated than the first generation, they explored how the 
timescale for spatial mixing of the two populations depends on 
the initial concentration of the second generation. They found that 
complete spatial mixing is expected only for dynamically evolved 
clusters which have lost at least 60 - 70% of their mass due to two- 
body relaxation (i.e. ignoring the rapid early loss of first-generation 
stars), irrespective of the initial concentration of the second gener¬ 
ation. In §5, we comment on the fraction of Galactic GCs expected 
to have evolved to that stage. 

In addition to having different spatial distributions, stars be¬ 
longing to different populations could also display distinct kine¬ 
matic properties. Bekki (2010, 2011) simulated the formation of 
second-generation stars from AGB ejecta within a cluster and 
showed that if the first generation formed with some small net 
angular momentum, the second generation can show considerable 
rotation and even adopt a flattened distribution initially. This is a 
simple consequence of conservation of angular momentum, along 
with dissipative processes driving the polluted material to higher 
densities towards the centre of the cluster. Such properties are rem¬ 
iniscent of the multiple stellar populations (with different ages and 
chemical properties) observed in nuclear star clusters, the youngest 
of which often form centrally concentrated stellar discs (e.g. Seth 
et al. 2006). Bekki (2010, 2011) also showed that the second- 
generation stars are expected to initially have a lower velocity dis¬ 
persion than the first-generation stars. 

Mastrobuono-Battisti & Perets (2013) were the first to study 
the long-term dynamical evolution of a GC with nested structures 
while also taking into account rotation and a flattened spatial dis¬ 
tribution for the more concentrated population. Their A-body sim¬ 
ulations, tailored to m Gen, trace the evolution of a relatively low- 
mass (10^ Mq) cold disc component embedded in a massive GC 
(2.5 X 10® Mq). They showed that these initial conditions can leave 
behind signatures even after 12 Gyr of relaxation-driven dynami¬ 
cal evolution'. For example, while the disc stars do become more 
isotropically distributed with time, they do not attain a completely 
relaxed spherical shape at the end of the simulations. Also, as these 
disc stars exchange angular momentum with the other stars, the 
whole cluster can become slightly flattened. Another prediction 
of this study is a lower velocity dispersion and a more radially 
anisotropic velocity distribution for the polluted (disc) stars com¬ 
pared to the pristine stars. We will address the uniqueness of these 
kinematic signatures in §4 when comparing the long-term evolution 
of the accretion and multiple generations scenarios. 

Before moving to our own A-body simulations, we note that 
the work presented here is also motivated by recent observational 
results. For example, from HST proper motion measurements in 47 
Tuc, Richer et al. (2013) found that the stars on the bluer side of the 
main sequence (presumably more helium-rich/polluted) exhibit the 
largest proper-motion anisotropy, with preferentially radial orbits, 
while the presumably pristine stars on the redder side of the main 
sequence display no velocity anisotropy (in the plane of the sky). 
The larger radial anisotropy of the polluted stars was interpreted 
by these authors as a possible effect of two-body relaxation, with 
the more centrally concentrated polluted population slowly diffus¬ 
ing outward. Kinematic differences were also observed between the 
stellar populations of 47 Tuc from radial velocity measurements. 
Kucinskas, Dobrovolskas & Bonifacio (2014) indeed reported that 


' Although note that the half-mass relaxation time of a cluster like ta Cen 
is comparable to a Hubble time or larger. 
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the line-of-sight velocity dispersion of their subsample of more en¬ 
riched stars is lower than that of their subsample of presumably 
pristine stars by about 1 km s“'. 

Based on radial velocities and Na abundance measurements in 
20 Galactic GCs, the study of Bellazzini et al. (2012) is the only one 
to date that has explored the possible connection between kinemat¬ 
ics and multiple chemical populations in a large sample of clusters. 
Their data generally did not reveal an obvious correlation between 
either velocity dispersion or systemic cluster rotation and Na abun¬ 
dance (a convenient tracer of self-enrichment). Three of their clus¬ 
ters (NGC 6388, NGC 6441, NGC 2808) show tentative evidence 
for a drop in the velocity dispersion with increasing Na abundance, 
but this will need to be verified with larger and more evenly dis¬ 
tributed samples. While it was reported that the Na-rich and Na- 
poor subsamples display similar rotational patterns in NGC 6441 
and NGC 6388, there is some indication that Na-poor stars dis¬ 
play a larger rotational amplitude than Na-rich stars in NGC 2808 
2.25 km s“' and 1.0 km s“', respectively^). Two other clusters 
(NGC 6171 and NGC 7078) show hints of Na-poor stars having a 
larger rotational amplitude by a few km s“*, albeit from small sam¬ 
ples. 


3 N-BODY SIMULATIONS 

3.1 Initial conditions for the pristine and polluted 
populations 

We studied the dynamical evolution of polluted and pristine stars 
in the context of the accretion (§3.1.1) and multiple generations 
(§3.1.2) scenarios by means of A-body simulations. More details 
about the aspects common to all of our simulations are summarised 
in the next subsection. Here we start by describing the initial con¬ 
ditions of the different populations and how these depend on the 
pollution scenario considered. 


3.1.1 Accretion scenario 


To capture the initial conditions of an accretion scenario like the 
one proposed by Bastian et al. (2013b), we set up an isochrone 
model (Henon 1959). The isochrone model is a convenient choice 
because the orbits in this potential are analytic, and we can easily 
flag as polluted the stars on an orbit crossing the core of the cluster. 
The potential of the isochrone model is given by 


Ofr) = - 


GM 

ro + (rl + r2)i/2 ’ 


( 1 ) 


where M is the cluster mass, r is the distance to the centre and tq 
is the scale radius (we shall refer to it as the core radius from now 
on), which is related to the half-mass radius of this model by r* = 
3.06 Tfl. The density profile, with a constant density core followed 
by a fall off (in 3D), is reminiscent of what is observed in nearby 
young massive clusters (Elson, Fall & Freeman 1987). 

To set up our systems, positions and velocities of the stars 
were sampled randomly from an isotropic isochrone model, and the 
masses of the particles were sampled from a Kroupa mass function 
(Kroupa 2001, see also §3.2), with no correlation between stellar 
mass and position. 


^ Bellazzini et al. actually quote values of Amt - 4.5 km s * and 2.0 
km s“*, but they define Arot as two times the mean rotational amplitude. 


We assumed that the enriched material - expelled by massive 
stars in the model proposed by Bastian et al. (2013b) - is located 
within the core radius of the cluster and we split the stars into two 
categories: stars that do not enter the core of the cluster are con¬ 
sidered as pristine stars, while stars on an orbit with a pericentre 
smaller than the core radius are assumed to accrete and are flagged 
as members of the polluted population. This yields a cluster con¬ 
taining about 45% of polluted stars at the start of the simulation, 
and these are by construction more centrally concentrated than the 
pristine stars. 

Note that both in Bastian et al. (2013b) and in the present 
work, the fraction of polluted stars is determined by the initial orbit 
distribution of low-mass stars. Early segregation (either primordial 
or dynamical) of the most massive stars is simply invoked as a way 
to get a source of polluted material in the central regions of the 
cluster. We do not set up our clusters with primordial mass segre¬ 
gation, but dense clusters that are not primordially mass segregated 
will quickly segregate by two-body interactions anyway ^. 

Not including primordial mass segregation could potentially 
affect the early dynamical evolution of the cluster, but we argue be¬ 
low that it does not change our results about the long-term evolu¬ 
tion. Early mass loss from centrally concentrated massive stars can 
lead to a stronger expansion than the same mass loss distributed 
throughout the body of the cluster. For tidally limited and strongly 
mass segregated clusters, Vesperini, McMillan & Portegies Zwart 
(2009) showed with N-body simulations that rapid dissolution may 
occur as a consequence of this early expansion and the associated 
larger flow of mass over the tidal boundary. However, they also 
showed that segregated clusters initially underfilling fheir Roche 
lobe survive the early expansion because much of this expansion 
occurs within the tidal radius and thus does not cause a significant 
loss of stars. In this case, the subsequent evolution, lifetime, and 
mass-loss rate of the cluster do not differ significanfly from those 
of an initially non-segregated cluster. The main notable difference 
is thaf long-lived initially segregated clusters will tend to have a less 
concentrated structure, delaying core collapse. Note that the early 
mass-loss driven expansion of a mass segregated cluster is not ho¬ 
mologous, as the massive (segregated) population loses more mass 
than the lower mass stars in the outer parts. The expansion is there¬ 
fore more important in the core and has much less severe effects in 
the outer parts. 

Because all our simulated clusters are initially tidally under¬ 
filling (as most GCs probably start) and because we are interested in 
long-term kinematic signatures that are more easily observed (and 
that more easily survive) in the outer parts of clusters, we conclude 
from the arguments above that our results would not be signifi¬ 
cantly affected if we had included primordial mass segregation. 

In the scenario proposed by Bastian et al. (2013b), low-mass 
stars may accrete a significant fraction of their final mass, which 
could modify the mass function of the polluted population. How¬ 
ever, we did not take this into account when building our N-body 
models for the accretion scenario, and instead assumed a standard 
mass function for both the polluted and the pristine population ini¬ 
tially. This is partly for the sake of simplicity, but also because cur¬ 
rently there are no clear predictions or robust observational con- 


^ A massive star of mass m sinks towards the centre of the cluster by dy¬ 
namical friction on a timescale f, = {m)lm X (Spitzer 1969), where 
is the relaxation time of the region of interest, and (m) is the mean stellar 
mass. The timescale fj can actually be shorter than a few Myr for dense 
clusters (Portegies Zwart et al. 1999). 
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straints on the mass function of the polluted population following 
this gas accretion phase. 

We only considered isochrone models with an isotropic veloc¬ 
ity distribution. Bastian et al. (2013b) showed that when consid¬ 
ering a radially anisotropic velocity distribution with the Osipkov- 
Merritt prescription (Osipkov 1979; Merritt 1985) and adopting an 
anisotropy radius equal to the half-mass radius, the fraction of stars 
that spend very little or no time in the core is in good agreement 
with the isotropic case. A radially anisotropic velocity distribution 
in the outer parts may represent more realistic initial conditions if 
clusters form via mergers of smaller clumps, but we will show in 
the next section that our simulated clusters rapidly develop radial 
anisotropy in the outer parts anyway. 

We computed various models with different amounts of ki¬ 
netic energy in rotation initially. This is motivated by studies re¬ 
vealing that rotation with a typical amplitude of a few km s“* is 
common among Galactic GCs (Cote et al. 1995; Anderson & King 
2003; van den Bosch et al. 2006; Lane et al. 2009,2010a,b; Bellazz- 
ini et al. 2012; Fabricius et al. 2014). Strong rotation has also been 
found in young and intermediate-age massive clusters (Henault- 
Brunet et al. 2012b; Mackey et al. 2013), while theoretical work 
suggests that cluster formation through violent relaxation in the 
presence of an external tidal field can lead to significant internal 
differential rotation (Vesperini et al. 2014). To quantify the initial 
amount of rotation of our models, we use the dimensionless spin 
parameter introduced by Peebles (1969): 


where J is the total angular momentum, E the total energy and M 
the mass of the system. For each scenario, we modelled clusters 
with 4 = 0, 0.091, and 0.129. In the specific case of an isochrone 
potential, these values of A correspond to fractions of respectively 
0, 10, and 20% of the total kinetic energy of the cluster in rotation. 
These are reasonable values considering the typical rotational am¬ 
plitudes and velocity dispersions of GCs (0 < Kot sin i /(XiD ^ 0.5; 
e.g. Meylan & Heggie 1997). For all our models with net rotation, 
the adopted direction of the rotation axis (the positive z-axis) is 
such that the cluster is corotating with its orbit around the centre of 
the galaxy. 

To consider clusters with different amounts of angular mo¬ 
mentum, we added rotation to the isochrone model by flipping 
the z-component of the angular momentum vector of a fraction 
of randomly selected stars having negative z-angular momentum 
within the cluster (which preserves virial equilibrium). The rotation 
curve initially has the same shape as the velocity dispersion pro¬ 
file, but as two-body relaxation proceeds it looks similar to what 
is obtained from self-consistent models including rotation (Varri 
& Bertin 2012), with a peak located near the half-mass radius of 
the cluster. We set up three systems (ACC0, ACCl, ACC2) with in¬ 
creasing values of A matching those used for the initial conditions 
of the multiple generations scenario (see Table 1 and §3.1.2). The 
stability of these systems was verified by checking that their total 
angular momentum, their half-mass radius, and their global virial 
ratio stayed constant for several crossing times after the start of the 
simulation (for test runs without stellar mass loss). 


3.7.2 Multiple generations scenario 

To emulate the initial conditions implied by the model of D’Ercole 
et al. (2008) and the hydrodynamical simulations of Bekki (2010, 


2011), we set up self-gravitating models with a flattened and rotat¬ 
ing stellar distribution embedded in a more extended spherical clus¬ 
ter. We preferred a purely dynamical approach over hydrodynami¬ 
cal simulations, partly to avoid the complications inherent to mod¬ 
elling the star formation process, but also because this approach is 
fast and gives us more control over the initial conditions. 

We used the code magalie (Boily, Kroupa & Penarrubia- 
Garrido 2001) included in the NEMO** Stellar Dynamics Tool¬ 
box (Teuben 1995, version 3.3.2). magalie can create multi- 
component (e.g. disc-fhalo) stellar systems in near equilibrium us¬ 
ing a method adapted from the BUILDGAL procedure of Hernquist 
(1993). When embedding a disc within a spherical component, the 
method starts with the exact spatial density profiles and then ap¬ 
proximates the phase space distribution function of all components 
(and hence the velocity field) using moments of the collisionless 
Boltzmann equation. For all our initial conditions generated with 
magalie (see below), we verified the stability of the constructed 
systems by checking that the angular momentum in both popula¬ 
tions, their half mass radii and the global virial ratio stayed con¬ 
stant for several crossing times and revolution periods of the disc 
after the start of the simulation (for test runs without stellar mass 
loss). 

We considered two distinct populations in our toy models, 
namely a polluted and a pristine one. For models with net angular 
momentum, we represented the pristine stars (first generation) by 
a spherically symmetric halo following a Hernquist profile^ (Hem- 
quist 1990) and the embedded polluted population (second gener¬ 
ation) by an exponential disc with a scale height equal to 20% of 
its scale length. Like the isochrone model, the Hernquist model has 
a power-law density profile with a fall off in the outer parts, 
which describes young massive clusters very well (Elson, Fall & 
Freeman 1987) and therefore represents a good choice for realistic 
initial conditions. 

We first explored how much angular momentum is expected 
to be lost from the first generation if ~ 90% of the initial cluster 
mass is rapidly removed from the outer parts. To do so, we set up 
systems in which the mass of the halo is 10 times larger than the 
mass of the disc (i.e. before mass loss). We assumed that half of the 
mass of the disc is made of ejecta from the first generation, and the 
other half comes from accretion onto the cluster of diluting (pris¬ 
tine) gas with zero net angular momentum. From conservation of 
angular momentum, we thus require that the total angular momen¬ 
tum of the halo is 20 times larger than that of the disc. The amount 
of rotation in the halo is controlled by flipping fhe z-componenf of 
fhe angular momentum vector of a fraction of randomly selected 
stars having negative z-angular momentum. A system satisfying 
the above constraints and having A = 0.091 is obtained by flip¬ 
ping the z-momentum of 65% of the halo stars having a negative 
z-momentum and by setting the ratio between the scale length of 
the Hernquist halo (flhaio) and the scale length of the exponential 
disc (T^di^c) fn n^haio/f^disc “ 9.37. 

Figure 1 shows the cumulative number and angular momen¬ 
tum distribution as a function of radius (in 3D) for the pristine and 
polluted populations of such a system. The vertical dashed line in 
each panel indicates the radius beyond which 90% of the cluster 
stars reside. If we assume that the stars beyond this radius are re¬ 
moved instantaneously, we can see that this would not only pref- 


http://bima.astro.umd.edu/nemo/ 

^ The density of the Hernquist model as a function of radius is given by 
= (r/a) a+r/a)3 ’ length. 
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Figure 1. Cumulative number and angular momentum distributions of pol¬ 
luted (blue curves) and pristine stars (red curves), as a function of radius, in 
an equilibrium “disc+halo” configuration with/i = 0.091, M^aio = lOMjjsc, 
and yz.haio = 20 /z.disc- This captures the configuration of the multiple gen¬ 
erations scenario before early violent mass loss. The vertical dashed lines 
indicate the radius beyond which 90% of the cluster stars reside. Top panel: 
Cumulative distribution of angular momentum (z-component) for the pol¬ 
luted (disc) and pristine (halo) stars as a function of radius, normalised to 
the total angular momentum in each component. Middle panel: Cumulative 
number distribution of the polluted and pristine stars as a function of radius, 
normalised to the total number of stars in each component. Bottom panel: 
Ratio of the number of polluted and pristine stars within radius r. 

erentially remove the pristine stars (as required in this scenario, 
leaving the number ratio of the two populations close to unity), but 
it would also remove a large fraction of the angular momentum of 
the pristine population. In this particular setup, the resulting cluster 
would have a larger total angular momentum in the second genera¬ 
tion than in the first generation by a factor of ~ 3. 

We also looked at similar configurations but with different val¬ 
ues of A and a larger initial mass ratio between the first and second 
generations (a mass ratio of 10 is a lower limit if the internal mass 
budget problem is to be avoided in scenarios with multiple genera¬ 
tions). For a fixed value of A, configurations with a larger mass ratio 
have a disc that is less compact with respect to the halo (i.e. smaller 
flhaio/f^disc)- In this case, the preferential loss of first-generation stars 
is less important when the outermost stars are removed. This can 
produce a system in which the total angular momentum of the first 
generation is still larger than that of the second generation, but it 
may also yield a fraction of second-generation stars that is too low 
(see Figure 2). As T is increased, ahaio/f^disc also gets smaller and 
it becomes increasingly difficult to remove a large fraction of pris¬ 
tine stars without also losing a large fraction of polluted stars. If 



% of total mass removed 


Figure 2. Top panel: Percentage of polluted stars in the cluster after in¬ 
stantaneous removal of all the stars beyond a certain radius (in the multiple 
generations scenario), as a function of the percentage of the total mass re¬ 
moved. Bottom panel: Ratio of total angular momentum in the remaining 
polluted and pristine stars, as a function of the percentage of the total mass 
removed from the cluster. In both panels, the different curves illustrate the 
effect of different initial conditions for the mass ratio between the two pop¬ 
ulations and the total amount of angular momentum (before mass removal). 


we take the simple dynamical considerations above at face value 
and if GCs can form with a large amount of angular momentum, 
a large initial mass ratio between the polluted and pristine popu¬ 
lations (Mhaio <: 20 Mjisc) rnay be difficult to reconcile with the 
observed fraction of polluted stars in GCs. 

In what follows, we focus on systems for which removing a 
large fraction of stars from the outer parts can leave the cluster 
with at least ~ 50% of polluted stars. In these cases, the total an¬ 
gular momentum in the remaining polluted stars is always larger 
than the total angular momentum in the remaining pristine stars 
(Figure 2). For simplicity, in all our models mimicking the mul¬ 
tiple generations scenario, we therefore assumed that the net an¬ 
gular momentum of the cluster is dominated by polluted stars and 
that the pristine population has no net angular momentum initially 
(•fz,pristine/-fz,polluted = 0). We also attributed an equal number of stars 
to each of the two populations. This captures the initial conditions 
after the early mass loss phase implied by this scenario. In the re¬ 
mainder of this paper, when we refer to mass loss, we are only 
concerned with the mass lost during the long-term evolution driven 
by two-body relaxation. 

The initial conditions of our different simulations are sum¬ 
marised in Table 1. For a subset of these (models MGENl, MGEN2, 
MGEN1 -nosev), we set up an exponential disc embedded in a spher¬ 
ical Hemquist halo. The spatial extent of the disc with respect to the 
halo is then dictated by the choice of A. One of these simulations 
(MGENl-nosev), with A = 0.091, is the one for which we did not 
include stellar evolution and which we followed until complete dis¬ 
solution (see §3.2). 
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Table 1. Summary of the initial conditions of our suite of W-body simulations. We also list the fraction of the initial mass and number of stars left after 10.75 Gyr 
for all the models with stellar evolution. The model without stellar evolution was followed until complete dissolution. More details on the parameters common 
to all models are presented in the main text. 


Model 

Density profile 

stellar 

evolution 

A 

^polluted /^pristine 

'^z.pristine/ polluted 

<3halo/^disc 

Mf/Mo 

Nf/No 

MGENl 

Exp. disc + Hemquist 

y 

0.091 

1 

0 

3.38 

0.48 

0.86 

MGENl-nosev 

Exp. disc + Hemquist 


0.091 

1 

0 

3.38 



MGEN2 

Exp. disc + Hemquist 

y 

0.129 

1 

0 

1.36 

0.50 

0.90 

Model 

Density profile 

Stellar 

evolution 

A 

^polluted /^pristine 

'^z.pristine/ polluted 

^pristine /^polluted 

Mf/Mo 

Nf/No 

MGENSa 

Hemquist + Hemquist 

y 

0 

1 


1.36 

0.46 

0.84 

MGENOb 

Hemquist + Hemquist 

y 

0 

1 


3.38 

0.45 

0.82 

MGENOc 

Hemquist + Hemquist 

y 

0 

1 


10.0 

0.42 

0.75 

Model 

Density profile 

Stellar 

evolution 

A 

^polluted /^pristine 

'^z.pristine/ ‘^z,polluted 

ro/rt 

A/f/il/o 

hlf/No 

ACCO 

Isochrone 

y 

0 

0.81 


0.33 

0.50 

0.90 

ACCl 

Isochrone 

y 

0.091 

0.81 

5.53 

0.33 

0.50 

0.90 

ACC2 

Isochrone 

y 

0.129 

0.81 

5.66 

0.33 

0.50 

0.90 


For comparison, we also considered another subset of initial 
conditions where the cluster has no net angular momentum and the 
second generation is represented by a spherically symmetric Flem- 
quist model embedded in a more extended Flemquist halo. Unlike 
the rotating models above, the spatial extent of the second genera¬ 
tion in these cases is not influenced by rotational support and does 
not depend on the properties of the first generation. Ultimately, it 
depends on the star formation process for this second generation. 
We thus decided to consider a few different concentrations for the 
second generation. In two cases (models MGENOa and MGENOb), the 
ratio of the radial scales of the two populations (apristine/apoiiuted) was 
taken to be the same as the flhaio/^disc values of our rotating models. 
We also considered a model in which the second generation is more 
concentrated (MGEN0C; see Table 1). 


3.1.3 Comparison of the initial conditions of the two scenarios 

As the polluted stars in the accretion model are assumed to be the 
ones crossing the core of the cluster, they are expected to be pref¬ 
erentially on radial orbits initially. In the presence of net angular 
momentum, this preference for radial orbits among the polluted 
stars also means that the mean rotational amplitude of these stars 
(at a given radius) will be lower than the mean rotational amplitude 
of the pristine stars. This is opposite to what is expected from the 
multiple generations scenario, where the angular momentum in the 
polluted population and the mean rotational velocity of the polluted 
stars is expected to be larger (Figure 3). 

Note that we have neglected the effect of gas accretion on the 
orbits of stars in our simulations. In the early disc accretion sce¬ 
nario, the enriched gas is released by massive stars near the centre 
of the cluster before being quickly swept up (within a few Myr). As 
these stars located near the centre have low orbital angular momen¬ 
tum compared to the rest of the cluster stars, the enriched material 
that they expel would also have low specific angular momentum. 
The specific angular momentum of a star would thus typically be 
reduced following accretion, and its orbit would shrink (through 
conservation of energy and angular momentum). If anything, tak¬ 
ing into account the effect of gas accretion would make the differ¬ 


ence in the mean rotational velocity (at a given radius) between the 
polluted and pristine populations even larger, enhancing the differ¬ 
ence already present between the initial conditions of the multiple 
generations and accretion scenarios. 

In Figure 3, we illustrate the similarities and differences in the 
initial conditions of the two scenarios by showing the number den¬ 
sity, velocity dispersion (z and x components), velocity anisotropy, 
and mean azimuthal velocity (around the z-axis) profiles at time 
t = 0 for the polluted and pristine populations of models MGENl 
and ACCl. Similar plots are shown in section A of the appendix to 
compare the initial conditions of models MGENOc and ACCO (Fig¬ 
ure Al) and models MGEN2 and ACC2 (Figure A2). 

The 3D number density profile is already very similar in both 
of the scenarios considered at the beginning of the simulation. The 
global density profile of all our models will evolve towards a King- 
like profile (King 1966) at late times, and any memory of the spe¬ 
cific model used to generate the initial conditions will gradually be 
erased as the density profile is reshaped by two-body relaxation and 
the tidal field. 

The initial conditions of both scenarios are also generally 
characterised by a dynamically cooler and more centrally concen¬ 
trated polluted population. Because the velocity dispersion of the 
cluster is higher towards the centre of the cluster, it may appear 
counter-intuitive that the more centrally concentrated population 
has a lower velocity dispersion. This is however a natural conse¬ 
quence of the more rapid fall off of the density distribution of this 
population (see e.g. van den Bosch et al. 1999). Another way to see 
this is that because the polluted population is preferentially located 
in the central regions, it contains fewer stars on wide orbits that 
cross the inner regions of the cluster at high velocity. 

For models MGENl and MGEN2, note however that the z compo¬ 
nent of the velocity dispersion is by definition initially smaller for 
the disc population than for the spherical (pristine) population. The 
only case where the velocity dispersion of the polluted population 
can be larger than that of the pristine population is when looking at 
the component of the velocity dispersion in the plane of the disc (cr^ 
or equivalently cty) for these two models, from about the half-mass 
radius and beyond. 

The anisotropy parameter is defined in the usual way as (3 = 
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MGENl; t=0 ACCl; t=0 



Figure 3. From top to bottom: number density, velocity dispersion in the z direction, velocity dispersion in the x direction, velocity anisotropy, and mean 
azimuthal velocity as a function of radius (in 3D) for polluted (blue) and pristine (red) stars. The initial conditions for model MGENl are shown in the left 
panels, and those for model ACCl are shown in the right panels. The half-mass radius is indicated by a grey dashed line. 


1 - ^, where cti and cr^ are the tangential and radial velocity dis¬ 
persions. The velocity distribution is thus radially anisotropic for 
j8 > 0. For all models with net rotation, the anisotropy profiles of 
the accretion and multiple generations scenarios are qualitatively 
similar. In both scenarios, the polluted population has a radially 
anisotropic velocity distribution, while the pristine population is 
either isotropic or even tangentially anisotropic in the inner regions 


for the early disc accretion scenario. Our accretion model with¬ 
out rotation (ACCO) displays the same behaviour, while both pop¬ 
ulations are (by construction) fully isotropic in the initial condi¬ 
tions of the multiple generations models without rotation (MGENQa, 
MGENOb, MGEN0C). 

The mean azimuthal velocity profile (also referred to as the 
“rotation curve” throughout this paper) is where the difference be¬ 
tween the initial conditions of the accretion and multiple gener- 
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t [Myr] 


Figure 4. Time evolution of the 1%, 5%, 10%, 25%, 50%, 75%, and 90% Lagrangian radii for the polluted (blue lines), pristine (red lines) and all stars (black 
lines) for all of our models with stellar evolution (see Table 1). 


ations scenarios is the most striking. Apart from models without 
rotation where the mean azimuthal velocity profile of both popula¬ 
tions is obviously flat and zero, all of our initial conditions imply 
significant differential rotation between the polluted and pristine 
stars. At a given radial distance from the centre of the cluster, the 
mean azimuthal velocity is always larger for the pristine stars in the 
accretion scenario, while it is always larger for the polluted stars 
in the multiple generations scenario. The latter follows from our 
choice of Tz.pristine/'/z,polluted = 0 for the initial conditions of the mul¬ 
tiple generations scenario. It would however remain true even with 
-/j,pristine/■/j,polluted = 1 because the polluted population would still be 
more centrally concentrated and would thus need to have a larger 
rotational velocity amplitude for its total angular momentum to be 
the same as the more extended pristine population. 


3.2 Description of the runs 


To follow the long-term dynamical evolution of the clusters, we 
used the code NB0DY6, a fourth-order Hermite integrator with Ah¬ 
mad & Cohen (1973) neighbour scheme (Makino & Aarseth 1992; 
Aarseth 1999, 2003) and accelerated force calculation on NVIDIA 
Graphical Processing Units (GPUs) (Nitadori & Aarseth 2012). 

All of our modelled clusters were placed on a circular orbit 
around the Galaxy, for which the gravitational potential was as¬ 
sumed to be a singular isothermal sphere. Unbound stars were re¬ 
moved from the simulation once they reached a distance of 2rj from 
the cluster centre, with rj (the Jacobi radius) given by 




( 3 ) 
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where Rq is the galactocentric distance, Vg is the circular velocity 
around the centre of the galaxy (we adopt Vg = 220 km s“')> G is 
the gravitational constant, and M is the mass of the cluster within 
radius rj. 

As we are interested in modelling mixing processes for which 
the evolution is driven by two-body relaxation (see Vesperini et al. 
2013), we did not include primordial binaries. This is a reason¬ 
able simplification because the post-core collapse expansion of the 
cluster and the accompanying evolution of the two-body relaxation 
timescale is expected to behave in the same way irrespective of the 
details of the source of energy driving the expansion (e.g. primor¬ 
dial binaries or mass loss from stellar evolution; Gieles, Heggie & 
Zhao 2011; Giersz & Heggie 2011). 

With one exception (see below), all of our models include stel¬ 
lar and binary evolution (and associated mass loss*) following the 
prescriptions of Hurley, Pols & Tout (2000) and Hurley, Tout & 
Pols (2002). All stars are given a metallicity of Z = 0.0035, but 
we ignore the effect of a spread in helium abundance among cluster 
stars on their evolution. For all our models with stellar evolution, 
we adopted a Kroupa (2001) initial mass function (IMF) between 
0.1 and 100 Mq, with a power-law slope of -1.3 below 0.5 Mq and 
-2.3 above. 

For the multiple generations scenario (§3.1.2), including stars 
up to 100 Mq may seem rather extreme given that we follow the 
evolution of the cluster after the second generation has formed and 
the cluster has reached virial equilibrium, in which case high-mass 
stars (> 8 Mq) from the first generation should have already ex¬ 
ploded as supernovae. By ignoring that first-generation massive 
stars have already exploded, we actually overestimate the early 
mass loss from supernovae by only ~ 10%. This may lead to 
an overestimate of the early expansion of the cluster by ~ 10% 
(e.g. Hills 1980), which is likely not very important when address¬ 
ing whether signatures of the initial conditions can survive for a 
timescale of the order of the relaxation time. 

For the model without stellar evolution, stellar masses were 
sampled from a double power-law IMF between 0.1 and 1.2 Mq (to 
mimic the present-day mass spectrum of GCs), with a break at 0.5 
Mq and power-law slopes of -1.3 and -2.0. 

All our models were computed with N* = 10* stars ini¬ 
tially. When including stellar evolution, the two-body relaxation 
timescale is tied to the timescale of mass loss from stellar evolu¬ 
tion, which fixes the physical timescale of the model. In these cases, 
the models with N’ = 10* stars can be scaled to represent a more 
massive cluster with a larger number of stars by demanding that the 
relaxation time is the same for the full-scale cluster with N stars as 
for the smaller-scale model. This is achieved when the radial scales 
of the two models are related in the following way (e.g. Heggie & 
Giersz 2008) 


_/A^y^V iog(yA^.) f' 
UJ \log(7A^)/ 


(4) 


where and represent the virial radius of the scaled and full- 
scale cluster, and y = 0.02 (e.g. Giersz, Heggie & Hurley 2008). 
The model with A* stars therefore has a larger radius than the full- 
scale cluster. We choose to scale our models to the properties of 
a typical massive GC like 47 Tuc, although we do not attempt to 
precisely reproduce the observed characteristics of this cluster. Our 
initial conditions are based on values from Giersz & Heggie (2011), 


who used Monte Carlo models to obtain a set of initial parameters 
providing a satisfactory match to the kinematic and photometric 
data of 47 Tuc after a Hubble time of dynamical evolution. At r = 0, 
our scaled models represent a cluster with A = 2 x 10* stars with 
a mean mass of (m) = 0.6377 Mq, a total mass of M = 1.275 x 
10* Mq, a virial radius of = 2.38 pc and a Jacobi radius of rj = 
86 pc, corresponding to a galactocentric radius of Rq = 3.3 kpc 
for a singular isothermal sphere potential and a circular velocity of 
Vg - 220 km s“'. The scaled parameters used were A, = 10*, total 
mass M, = 6.377 x 10"* Mg, /•„ , = 5.2 pc and Rq ,, = 47.98 kpc. 

We evolved these clusters up to an age of 10.75 Gyr, which 
is close to the latest age determination of 47 Tuc from the proper¬ 
ties of its white-dwarf cooling sequence (Hansen et al. 2013). At 
the end of the simulations, the modelled clusters are all reason¬ 
ably consistent with the present-day characteristics of 47 Tuc (see 
Giersz & Heggie 2011): half-mass radii between 5 and 7 pc, central 
line-of-sight velocity dispersions between 10 and 13 km s“', Jacobi 
radii between 64 and 69 pc. Table 1 lists the fraction of the initial 
mass (Mf/Mo) and number of stars (Af/Ao) left at the end of the 
simulation for all these models with stellar evolution. 

For the model without stellar evolution, a cluster with 
10* stars, a virial radius of 2 pc and a mean mass (m) = 0.3458 Mg 
was placed on a circular orbit with Rq = 10 kpc. Because the 
physical timescale of this model is not set by the stellar evolution 
timescale, we are free to scale both the mass and radius of the sys¬ 
tem to represent clusters with different relaxation times and galac¬ 
tocentric radii We ran this model until complete dissolution of the 
cluster, which allows us to explore the late stages of the dynamical 
evolution of multiple populations without being tied to a specific 
scaling. 

The initial half-mass relaxation time is 0.9 Gyr for all our sim¬ 
ulations with stellar evolution and 0.3 Gyr for our model without 
stellar evolution, using the equation from Spitzer (1987), which as¬ 
sumes an equal-mass cluster and spherical symmetry. One should 
however be cautious when directly comparing the age of the sys¬ 
tems to this initial half-mass relaxation time, first because the half¬ 
mass relaxation time is not well defined for systems with a disc em¬ 
bedded in a spherical halo. We should expect two-body relaxation 
to be more efficient in systems in which stars are moving coherently 
(e.g. the disc component of our multiple generations systems), but 
this is not captured by the equation used. Our clusters also have a 
stellar mass spectrum and the mass function changes with time as 
a result of stellar and dynamical evolution. This will have an effect 
on the two-body relaxation time at any given point of the evolution. 


4 LONG-TERM KINEMATIC IMPRINTS 
4.1 Phase-space mixing 

We now turn to the long-term dynamical evolution of clusters 
from the initial conditions described in the previous section. Af¬ 
ter 10.75 Gyr, our 47 Tuc-like models with stellar evolution have 
lost between 50% and 58% of their initial mass (~ 10-25% of their 
stars) through a combination of stellar mass loss and tidal evapora¬ 
tion. At this stage, preferential loss of pristine stars has increased 
the ratio of polluted to pristine stars to Apoii„ted/Npristine ~ 0.85 for 
the accretion models and 1.03 < ApoUuted/Apnstine < 1-42 for the 
multiple generations models. 


* The mass lost due to stellar evolution is assumed to be instantly removed 
from the cluster. 


^ Note that the ratio of the half-mass radius to the Jacobi radius (ri,/rj) is 
fixed. 
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Figure 5. Evolution of the 1%, 5%, 10%, 25%, 50%, 75%, and 90% La- 
grangian radii for the polluted (blue lines), pristine (red lines) and all stars 
(black lines) as a function of the fraction of the initial mass left in the clus¬ 
ter for model MGENl-nosev (see Table 1). The radius scale of this scalable 
model is expressed in W-body units. 


Figure 4 shows the time evolution of the Lagrangian radii of 
polluted and pristine stars (along with those of the whole cluster) 
for all the models with stellar evolution. Overlap of all Lagrangian 
radii for the two populations would indicate complete spatial mix¬ 
ing, but for none of the models in Figure 4 are the polluted and pris¬ 
tine stars fully mixed at the end of the simulation. Only in model 
MGENOa are both populations almost mixed, but this is the model 
in which the polluted population started the least concentrated with 
respect to the pristine population. Interestingly, the effect of the 
gravogyro instability in accelerating core collapse (e.g. Ernst et al. 
2007) is seen when comparing models ACCO, ACCl and ACC2, for 
which the initial conditions only differ in their total amount of an¬ 
gular momentum. 

Vesperini et al. (2013) showed that both internal two-body re¬ 
laxation and mass loss drive the evolution of the mixing of multiple 
populations. The relaxation timescale is longer in the outer regions 
of the cluster, making two-body relaxation (and mixing) less effi¬ 
cient at these larger distances from the centre. This is amplified by 
the expansion of the cluster which slows down the rate of dynam¬ 
ical ageing with time. Mass loss acts to slow the growth of the re¬ 
laxation time, which eventually starts to decrease as the remaining 
mass and the radius of the cluster decrease. The Jacobi radius also 
decreases as the cluster mass goes down, and the less mixed outer 
layers are gradually stripped away, which accelerates the evolution 
towards complete spatial mixing. 

Regardless of initial differences in the concentration of the 
polluted population, Vesperini et al. (2013) found that clusters al¬ 
ways approach a state of complete mixing after losing 60 - 70% 
of their mass. This is also what we find in model MGENl-nosev, 
which was followed until complete dissolution. Figure 5 shows the 
evolution of the Lagrangian radii of the polluted and pristine stars 
for this model as a function of the fraction of the initial mass left 
(MIMq). All Lagrangian radii for the two populations overlap when 
M/Mo ~ 0.35. 

In the present work, we are mainly interested in the kinematic 
properties of multiple populations as a way to distinguish between 
the proposed scenarios. It is therefore important to verify if these 
kinematic imprints are erased on the same timescale as full spatial 



Figure 6. Distribution of energy and z-angular momentum for the polluted 
(blue) and pristine (red) stars of model ACCl at r = 0 (top panel) and at the 
end of the simulation (bottom panel). The solid lines represent isodensity 
contours in the |£| vs. J. plane for each population. 


mixing occurs. In their study of the long-term dynamical evolution 
of globular clusters with two stellar populations (one of them be¬ 
ing more concentrated initially), Decressin, Baumgardt & Kroupa 
(2008) found that the loss of information on the stellar orbital angu¬ 
lar momentum occurs on the same timescale as spatial homogeni¬ 
sation. To check this, we looked at the evolution of the two popula¬ 
tions of our models in energy and angular momentum space. 

Figure 6 shows the distribution of energy and z-angular mo¬ 
mentum for the two populations of model ACCl at the beginning 
and at the end of the simulation. The energy of a star is defined as 
E = tfnp + m(l>c, where (pc - the potential energy due to all other 
cluster stars - excludes the contribution of the tidal field. One can 
see the obvious differences in the initial conditions of the two pop¬ 
ulations for the accretion scenario. Polluted stars typically having 
a lower energy (i.e. more negative - they are more tightly bound to 
the cluster due to their smaller distance from the centre and their 
low velocity). The net angular momentum of this model and the 
larger angular momentum of the pristine population are also appar¬ 
ent from the distribution of J-. After 10.75 Gyr, differences are still 
present both in the energy and angular momentum distributions of 
the two populations. 

Figure 7 shows the same quantities but for model MGENl. 
Again the polluted stars typically have a lower energy. The angular 
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Figure 7. Distribution of energy and z-angular momentum for the polluted 
(blue) and pristine (red) stars of model MGENl at t = 0 (top panel) and at the 
end of the simulation (bottom panel). The solid lines represent isodensity 
contours in the |£| vs. J. plane for each population. 


momentum distribution of the pristine stars is initially symmetric 
around = 0, while the polluted stars have an angular momen¬ 
tum distrihution skewed towards positive values of J, for this mul¬ 
tiple generations model with net angular momentum. After 10.75 
Gyr, the two populations are again not fully mixed in energy and 
angular momentum space. Some exchange of angular momentum 
has occurred between the polluted and pristine stars, as the final 
distribution of pristine stars is slightly skewed towards positive 
values*. 

As pointed out by Mastrobuono-Battisti & Perets (2013), par¬ 
tial relaxation and the exchange of angular momentum between the 
polluted and pristine populations can have consequences on the 
spatial morphology of the two populations (and the cluster as a 
whole) when the polluted population starts as a disc/flattened struc¬ 
ture. We illustrate this in Figure 8 hy showing the cumulative frac¬ 
tion of stars as a function of cos i (where i is defined as the position 
angle of the star with respect to the positive z-axis) for the polluted, 
pristine, and all stars of model MGENl at the beginning and at the 
end of the simulation. The polluted population has not relaxed to a 


° Although we note that tidal effects can also induce some mild rotation. 
We will come back to this briefly in §4.4 


MGENl 

t=0.0 Myr / 



- pristine 

- polluted 

- all 



Figure 8. Cumulative fraction of stars as a function of cos i (where i is the 
position angle of the star with respect to the positive z-axis) for the polluted 
(blue), pristine (red) and all stars (black) of model MGENl at f = 0 (top 
panel) and at the end of the simulation (bottom panel). Stars in the plane of 
the disc (the x — y plane) have cos i = 0. The dashed grey line represents an 
isotropic spatial distribution. 


completely spherical shape (the isotropic distribution is represented 
by the dashed grey line in Figure 8) hy the end of the simulation. 
Slight flattening of the pristine population (which started isotropic) 
and of the cluster as a whole also follows from the exchange of an¬ 
gular momentum between the two populations, although the pris¬ 
tine population is obviously still less flattened that the polluted one. 
Studying the ellipticity profiles of multiple populations in globular 
clusters may reveal important clues about their formation. To es¬ 
tablish whether this would allow to distinguish between different 
scenarios, one should however also consider accretion models in 
which the cluster is initially flattened, but this is beyond the scope 
of the present work. 

Like the two examples shown in Figures 6 and 7, none of 
our 47 Tuc-like models are fully mixed in energy and angular mo¬ 
mentum space after 10.75 Gyr (see figures in section B of the ap¬ 
pendix). To estimate when full mixing occurs, we turn to model 
MGENl-nosev again. Figure 9 shows the distrihution of angular 
momentum and energy of polluted and pristine stars in this model 
at three different moments of the evolution of the cluster: when 
M/Mo = LO, 0.5, and 0.35. After the cluster has lost 50% of 
its initial mass, the distrihution of angular momentum of the pol¬ 
luted stars is still skewed towards larger positive values of J, com¬ 
pared to that of the pristine stars, similar to what was found for 
model MGENl. The energy distribution of the polluted stars is also 
skewed towards lower (more negative) values compared to the pris¬ 
tine stars. When M/Mo = 0.35, the shapes of the liil and distribu¬ 
tions are essentially the same for polluted and pristine stars (up to a 
scaling factor due to the different number of stars in the two popu¬ 
lations). This suggests that kinematic differences befween polluted 
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Figure 9. Distribution of energy and z-angular momentum for the polluted 
(blue) and pristine (red) stars of model MGENl-nosev at t = 0 (top panel), 
when M/Mo = 0.50 (middle panel) and when M/Mq = 0.35 (bottom panel). 
The solid lines represent isodensity contours in the |£| vs. plane for each 
population. 


and pristine stars are erased on the same timescale as full spatial 
mixing occurs. 

In the next subsections, we look in more details into the spe¬ 
cific kinematic signatures that are expected from each scenario 
when a cluster is not fully mixed. 


4.2 Velocity dispersion 

We now briefly discuss the evolution of the velocity dispersion pro¬ 
file of the two populations for our models with stellar evolution. 
Figure 10 shows the time evolution of the velocity dispersion (the 
z-component) profile of the polluted and pristine stars for models 
MGENl and ACCl. 

The variation of the two-body relaxation time as a function 
of radius and time is clearly seen in Figure 10. After only about 
1 Gyr, the central z-velocity dispersion is the same for the polluted 
and pristine stars in both models. Nearly 10 Gyr later, the region 
within which the velocity dispersion profile of the two populations 
is identical extends further out, but differences of the order of 1 or 2 
km s“' persist in the outer parts of the cluster where the relaxation 
time is longer. The region where these differences persist appears 
to extend further in for model MGENl, but when also considering 
simulations with other values of A and different concentrations of 
the polluted population initially (see Figures Cl and C2 in the Ap¬ 
pendix), this is not always apparent. It does not appear to be an un¬ 
ambiguous signature that would allow to distinguish between the 
scenarios. 

Note that for a dynamically young system, differences in the 
velocity dispersion of the polluted and pristine stars may persist 
even in the inner regions of the cluster. This could explain why the 
core velocity dispersions of the first and second generation stars 
are still different after 12 Gyr in the simulations of Mastrobuono- 
Battisti & Perets (2013), as their disc-Ehalo model is tailored to 
a) Cen. 

As discussed by Bekki (2010, 2011) and Mastrobuono-Battisti 
& Perets (2013), a lower velocity dispersion for the polluted pop¬ 
ulation is one kinematic signature that can be left in a scenario 
with multiple generations, especially when looking at the compo¬ 
nent perpendicular to the initial plane of the disc. This is however 
clearly not unique to such a model, as we find the same signature 
in the accretion scenario, regardless of the initial value of A. There¬ 
fore, the lower velocity dispersion reported for the enriched stars in 
47 Tuc (Richer et al. 2013; Kucinskas, Dobrovolskas & Bonifacio 
2014) cannot be used to favour or discard one of the two scenarios 
considered here. 

As mentioned in §3.1.3, the initial x (ory) velocity dispersion 
(i.e. the component in the rotation plane of the disc) of the polluted 
population can be larger than that of the pristine population in a 
multiple generations model with net angular momentum. Figure 11 
shows the time evolution of the velocity dispersion (the x compo¬ 
nent) profile of the polluted and pristine stars for models MGENl 
and ACCl (see Figures C3 and C4 in the Appendix for cases with 
different values of A initially). 

At f = 10.75 Gyr, there is almost no trace of the initially larger 
X velocity dispersion of the polluted population in model MGENl, 
as both the polluted and pristine populations have essentially the 
same x dispersion profile. In model MGEN2 (Figure C4), a larger 
X velocity dispersion for the polluted population is however still 
apparent at the end of the simulation in the outer parts of the cluster, 
but the difference is very small (< 1 km s“*). 

For the accretion models and the multiple generation models 
without a disc (i.e. no net angular momentum), which are all spheri¬ 
cally symmetric, the x and y components of the velocity dispersion 
are of the same magnitude as the z component, so the arguments 
above about the z dispersion still hold for these systems. The veloc¬ 
ity dispersion (any component) is always smaller for the polluted 
(and more concentrated) population in these cases. 
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ACCl 

t=1386 Myr 




Figure 10. Velocity dispersion (z-component) as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGENl (left panels) and 
model ACCl at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines. 


4.3 Velocity anisotropy 

Motivated by the finding of a radially anisotropic polluted popula¬ 
tion in 47 Tuc (Richer et al. 2013), we now look at the evolution of 
the velocity anisotropy profile of the polluted and pristine stars in 
the accretion and multiple generations scenarios. Figure 12 shows 
the time evolution of this profile for each of the two populations in 
models MGENl and ACCl. 

In both models, the polluted population stays more radially 
anisotropic than the pristine population in the outer parts of the 
cluster. The anisotropy profile of the polluted stars is characterised 
by an inner isotropic core followed by radial anisotropy that in¬ 
creases with radius. In the late stages, the radial anisotropy peaks 
at an intermediate distance from the centre and then decreases out¬ 


wards, as expected when the effect of tides becomes important 
and stars on radial orbits are preferentially lost (e.g. Baumgardt 
& Makino 2003; Fukushige & Heggie 2000) and gain angular mo¬ 
mentum due to the tides (Oh & Lin 1992). After 10.75 Gyr, the 
pristine population is also isotropic in the inner parts (over a larger 
radial extent than the polluted population) and shows some very 
mild radial anisotropy at larger radii, followed by a decrease of 
and isotropy or mild tangential anisotropy in the outermost regions. 

These conclusions do not change in models with different val¬ 
ues of A. Even when 4 = 0 and both populations of a model 
with multiple generations start fully isotropic across the whole 
cluster, the more rapidly expanding (and initially more centrally 
concentrated) polluted population rapidly develops stronger radial 
anisotropy in the outer parts compared to the pristine population. 
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Figure 11. Velocity dispersion (x-component) as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGENl (left panels) and 
model ACCl at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines. 


The kinematic signature at the end of the simulation in this case is 
remarkably similar to all of our other models (Figure 13). A more 
radially anisotropic polluted population is therefore not a unique 
signature of either scenarios, but another natural consequence of 
starting with a more centrally concentrated polluted population. 


4.4 Rotation 

We now look at the long-term evolution of the bulk rotation of 
pristine and polluted stars, a promising diagnostic given the very 
distinct initial conditions of the accretion and multiple generations 
scenarios for this kinematic property (Figure 3). To illustrate the 
differences between the two scenarios as a function of time, Fig¬ 


ure 14 shows the evolution of the mean azimuthal velocity profile 
for the two populations in models MGENl and ACCl. 

As the clusters evolve due to two-body relaxation and lose 
mass, angular momentum is transported outwards and the total an¬ 
gular momentum decreases. The peak of the rotation curve (located 
near the half-mass radius) also moves outward and its amplitude de¬ 
creases. In all models with net rotation, as mixing progresses, the 
mean azimuthal velocity of the pristine stars eventually matches 
that of the polluted stars in the inner regions. The rotational am¬ 
plitude in these inner regions is very low after 10.75 Gyr, but a 
small velocity gradient (of the order of 1 km s“' pc“') is still present 
across the core of the cluster, even in post-core collapse clusters (for 
example model MGENl, see Figure 4). This suggests that observa¬ 
tional evidence for a velocity gradient/rotation in the inner parts of 
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Figure 12. Velocity anisotropy (jS) as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGENl (left panels) and model ACCl 
at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines. 


a cluster cannot be used to argue that the cluster is in the pre-core 
collapse phase of its evolution (c.f. Fabricius et al. 2014). 

In all our models with net rotation, the differential rotation of 
the pristine and polluted stars is maintained in the outer parts of 
the cluster after 10.75 Gyr of dynamical evolution. For the mul¬ 
tiple generations scenario, the polluted stars rotate faster than the 
pristine stars, while the opposite is true for the accretion scenario. 
While the rotation curves of the two populations become indistin¬ 
guishable in the inner regions of the clusters, differences of the or¬ 
der of 1 - 2 km s“* or more persist in the outer parts^. A larger 


^ As discussed in §3.1.1, taking into account the effect of gas accretion on 


initial value of A yields larger differences in the rotation curves of 
pristine and polluted stars (see Figure D2 in the appendix). 

In models MGENl and MGEN2 (Figures 14 and D2), the ex¬ 
change of angular momentum from the polluted to the pristine pop¬ 
ulation is seen through the subtle growth of the rotational amplitude 
of the pristine population (which has zero angular momentum ini¬ 
tially) in the inner parts of the clusters. The rotation curve of the 
pristine stars also rises very slightly in the outer parts of these clus¬ 
ters as they evolve, showing a gain of weak prograde rotation at 
the level of ~ 1 km s“'. This effect is also seen in our models with 


individual stellai' orbits would likely enhance the amplitude of the differen¬ 
tial rotation between pristine and polluted stars in the accretion scenario. 
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r [pc] r [pc] 

Figure 13. Velocity anisotropy (j8) as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGENSc (left panels) and model ACCS 
at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines. 


no initial rotation (Figure Dl), but the effect remains small and is 
always confined to regions beyond the tidal radius of the cluster. 
Tidal effects are expected to cause a preferential loss of stars on 
prograde orbits (Henon 1970; Keenan & Innanen 1975; Oh & Lin 
1992; Fukushige & Fleggie 2000), so we are most likely witnessing 
the stripping of stars on prograde orbits when seeing an increas¬ 
ing prograde rotation outside of the tidal radius. In any case, for 
our models with initial rotation, such a tidally induced rotational 
signature never compromises the clear differential rotation signal 
between the polluted and pristine stars. 


5 PROSPECT OF DETECTING THE IDENTIFIED 
SIGNATURES 

Based on the results of the previous section, comparing the veloc¬ 
ity dispersion of the polluted and pristine populations may offer a 
way to distinguish between the scenarios kinematically, although 
the expected signal is weak. A larger velocity dispersion for the 
polluted population in the plane of rotation can be explained by a 
multiple generations model with net angular momentum, but not 
by the accretion model nor any model without net rotation initially. 
Observationally, this component of the velocity dispersion is best 
probed with radial velocities if the rotation axis is perpendicular 
to the line of sight, but best probed with proper-motion data if the 
rotation axis is along the line of sight. 
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Figure 14. Mean azimuthal velocity as a function of radius (in 3D) for the polluted (blue) and pristine (red) stars of model MGENl (left panels) and model 
ACCl at different times in the evolution of the cluster. The half-mass radius is indicated by grey dashed lines. The inset plot in the bottom right panel shows a 
zoomed view of a portion of the rotation curve for the final snapshot of model ACCl. 
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M{t)/Mo 

Figure 15. Amplitude of the rotation curve (i.e. peak height of the mean 
azimuthal velocity profile) for the polluted (blue lines), pristine (red lines) 
and all stars (black lines) as function of the fraction of the initial mass left 
in the cluster for model MGENl-nosev. 


If a larger velocity dispersion is measured for the polluted pop¬ 
ulation in a given dataset, the multiple generations model would be 
favoured and the same data should reveal a larger rotational ampli¬ 
tude for this polluted population. On the other hand, if a lower ve¬ 
locity dispersion is measured for the polluted population, it could 
mean that the accretion model prevails but could also be because 
the inclination of the rotation axis is not favourable enough to probe 
the velocity dispersion in the plane of rotation. In that sense, the ve¬ 
locity dispersion is perhaps not the most powerful probe of the for¬ 
mation scenario. On the other hand, it is clear that the differential 
rotation of polluted and pristine stars represents a unique kinematic 
signature that can help to constrain scenarios for the formation of 
multiple populations in GCs. 

The expected strength of the signal in a cluster like 47 Tuc 
with a moderate amount of angular momentum initially is a few 
km s“* or less. This could be observable with precise ground-based 
radial velocity measurements from large high-resolution spectro¬ 
scopic samples of RGB stars (bearing in mind inclination effects). 
As discussed in §2, tantalising results in that respect have been ob¬ 
tained by Bellazzini et al. (2012), but they will need to be confirmed 
with larger datasets or at least with a more in-depth statistical analy¬ 
sis. Interestingly, in the three clusters for which these authors report 
different rotational amplitudes for the Na-rich and Na-poor sam¬ 
ples, the pristine (Na-poor) stars always display a larger rotational 
amplitude, as expected in the accretion scenario. 

The Gaia satellite will also allow to obtain proper motions 
with a precision better than 2 km s“' for RGB stars in the out¬ 
skirts of GCs, which could be sufficient to see a difference in the 
rotation of polluted and pristine stars in some clusters. Gaia will 
not be helpful for proper motions in the dense cores of GCs, but 
for the application discussed here, this is not a problem because the 
signal is expected to be stronger in the outer parts of clusters, where 
the two-body relaxation timescale is longer. Proper motions from 
HST (e.g. Richer et al. 2013) may also be useful to uncover the 
expected signal. In this case, the analysis would be made easier by 
the differential nature of the signal, which eliminates the need for 
determining an absolute reference frame and avoids the associated 
complications, especially when the field-of-view of the observa¬ 
tions is small (e.g. Anderson & King 2003). 

To estimate the fraction of the Galactic population of GCs in 
which we can hope to detect differential rotation of the pristine and 
polluted stars, we must identify the clusters for which the different 
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Figure 16. Distribution of stellar mass and galactocentric radius of Galac¬ 
tic GCs from the 2010 version of the Harris catalogue. Clusters above the 
dashed line are presumably in the expansion-dominated phase of their evo¬ 
lution and not fully mixed, while clusters below the line are expected to 
be in the evaporation-dominated phase. 47 Tuc is represented by a green 
square, and the fully spatially mixed cluster NGC 6362 by a red cross. 


populations are not expected to be fully mixed. As we argued pre¬ 
viously, these are the clusters which have lost less than ~ 60 - 70% 
of their initial mass (during the long-term evolution driven by two- 
body relaxation). Figure 15 shows that the amplitude of the rota¬ 
tion curve is the same for the polluted and pristine stars of model 
MGENl-nosev when the cluster has lost 70% of its mass or more, 
supporting the idea that the differential rotation signature would be 
erased beyond that stage. 

To identify the parameters of the clusters which are more 
likely to have lost a large fraction of their mass, we follow Gieles, 
Heggie & Zhao (2011), who presented a simple description of the 
life cycle of initially compact clusters in a tidal field. They showed 
that the half-mass radius of a cluster increases during (roughly) 
the first half of its evolution (the so-called ‘expansion dominated’ 
phase) and decreases during the second half (the ‘evaporation dom¬ 
inated’ phase). Clusters in the evaporation-dominated phase have 
lost about 50% of their mass or more. From their evolution equa¬ 
tions and adopting a few simplifying assumptions (e.g. circular 
cluster orbits in a constant Milky Way potential), Gieles, Heggie 
& Zhao (2011) showed that clusters satisfying the relation 


M > 10= Mq 


/ 4 kpc \ 


(5) 


should (to first order) be in the expansion-dominated phase, in 
which case we do not expect them to be fully mixed. In Fig¬ 
ure 16, we show clusters in the M - Rq plane using data from 
the 2010 version of the Harris catalogue*** (Harris 1996), where we 
adopted a mass-to-light ratio of 2 to convert luminosities to masses 
(McLaughlin & van der Marel 2005). There is a large sample of 
Galactic GCs (roughly half of the data points) with large enough 
masses and galactocentric radii to belong to the group of expansion- 
dominated clusters. This is where the kinematic signature identified 
in the present work should be looked for. 47 Tuc is found in this re- 


*** http://www.physics.mcmaster.ca/Globular.html 
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gion of the M - Rq plane, which is consistent with the kinematic 
differences reported between its pristine and polluted stars (Richer 
et al. 2013; Kucinskas, Dohrovolskas & Bonifacio 2014). On the 
other hand, the lower-mass cluster NGC 6362, the first GC found 
to he fully spatially mixed (Dalessandro et al. 2014), is interest¬ 
ingly very close to the boundary between the expansion-dominated 
and evaporation-dominated regimes. 


6 CONCLUSIONS 

We focused on identifying a unique long-term kinematic imprint 
that would allow to distinguish between different scenarios for how 
enriched material makes its way into stars in the context of mul¬ 
tiple stellar populations in GCs. To do so, we presented a suite of 
Al-body simulations of the long-term dynamical evolution of glob¬ 
ular clusters, with initial conditions chosen to capture the distinct 
kinematic properties of two main pollution scenarios. 

The first scenario (the multiple generations scenario) is one 
in which gas collects in a cooling flow into the core of the cluster, 
where a new generation of stars is eventually formed. In the other 
scenario (the accretion scenario), a fraction of the stars accrete ma¬ 
terial when passing through the core of the cluster, with no need 
for multiple star-formation events. While in the current paradigm 
these scenarios are naturally associated to the models put forward 
by D’Ercole et al. (2008) and Bastian et al. (2013b) and the specific 
source of polluting material favoured by these authors (AGB stars 
and massive interacting binaries, respectively), our study does not 
concern the chemistry and our results would apply to scenarios that 
proceed in the same way dynamically, regardless of the source of 
pollution. 

We showed that the two scenarios imply different initial con¬ 
ditions for the kinematics of the various populations. The most 
striking difference befween the initial conditions of the accretion 
and multiple generations scenarios arises in the presence of cluster 
rotation. As a result of dissipative processes driving the polluted 
material to higher densities towards the centre of the cluster and 
conservation of angular momentum, an initially larger rotational 
amplitude for the polluted stars compared to the pristine stars is 
expected from the multiple generations scenario. In the accretion 
scenario, the polluted stars (the ones crossing the core) are pref¬ 
erentially on radial orbits and their initial rotational amplitude is 
expected to be smaller than that of the pristine stars. 

We showed that initial differences in the kinematics of pris¬ 
tine and polluted stars can survive for a Hubble time of relaxation- 
driven dynamical evolution in a cluster with properties similar to 
those of 47 Tuc, especially in the outer parts of the cluster where the 
two-body relaxation timescale is longer. In particular, some differ¬ 
ential rotation of the pristine and polluted stars is preserved, which 
means that this unique signature can potentially be used to dis¬ 
tinguish between the pollution scenarios considered. The expected 
strength of the differential rotation signal is typically a few km s“' 
in the outer parts (the exact value depends on the specific scenario, 
the initial amount of angular momentum of the cluster, the radius 
at which this is measured, and on the degree of mixing - i.e. the 
dynamical age of the cluster). This is challenging to measure but 
within the reach of current and future surveys of the internal kine¬ 
matics of GCs (either through radial velocities or proper motions). 

Initial differences in the velocity dispersion and velocity 
anisotropy profiles of polluted and pristine stars are also expected 
to survive for a Hubble time in a 47 Tuc-like cluster. These kine¬ 
matic imprints are however remarkably similar in the accretion and 


multiple generations scenarios. A lower velocity dispersion and a 
more radially anisotropic velocity distribution for the polluted stars 
(common signatures of both scenarios) can therefore not be used 
to distinguish the two scenarios. That said, the multiple genera¬ 
tions scenario could be favoured if a larger velocity dispersion is 
measured for the polluted population. This would stem from the 
initially larger velocity dispersion of this population in the plane of 
rotation. 

We argued that full kinematic mixing of the two populations 
happens on the same timescale as full spatial mixing. It is achieved 
when ~ 60 - 70% of the initial mass of the cluster has been 
lost. We identified a large subsample of Galactic GCs having large 
enough masses and galactocentric radii to still be in the expansion- 
dominated phase of their evolution (Gieles, Heggie & Zhao 2011). 
Clusters in this phase are still expanding within their tidal bound¬ 
ary and have not suffered substantial tidal evaporation. They are 
thus good candidates for observing the signatures of incomplete 
spatial and kinematic mixing, in particular the differential rotation 
of pristine and polluted stars. 

We should mention that the present work does not encompass 
all suggested frameworks for the formation of multiple populations. 
For example. Maxwell et al. (2014) suggested that GCs formed in 
the inner regions of dwarf galaxies can accrete gas on each passage 
through the centre of these galaxies. The accreted matter, eventu¬ 
ally forming a second generation within the cluster, would in this 
case be a mix of pristine material from the intergalactic medium 
and surrounding dwarf, plus moderate velocity AGB ejecta retained 
by the potential well of the galaxy. These GCs would be pushed out 
onto larger orbits as star formation progresses, and may be stripped 
during mergers. While it avoids some of the problems of the other 
models, this scenario cannot explain the presence of multiple stellar 
populations in the Galactic GCs thought to form in situ (e.g. Forbes 
& Bridges 2010; Leaman, VandenBerg & Mendel 2013). It would 
however be interesting to see if this scenario can leave a charac¬ 
teristic kinematic signature. Similarly, we have not investigated the 
implications of another scenario in which second-generation stars 
are proposed to form in the decretion discs of fast rotating massive 
stars (Krause et al. 2013). 

Finally, there is definitely scope to improve the toy models that 
we have set up to study the kinematics of multiple populations in 
different scenarios. For example, the effect of gas accretion on the 
orbits of stars could be included in a self-consistent way. It would 
also be interesting to consider the effect of the different initial he¬ 
lium content of the polluted and pristine stars on their evolution, 
and indirectly on the dynamical evolution of the cluster. 
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APPENDIX A: SUPPLEMENTARY FIGURES (INITIAL 
CONDITIONS) 

This section presents two additional figures to complement §3.1 
and Figure 3 and illustrate initial conditions for models without 
net angular momentum initially (MGENOc and ACC0; Figure Al), 
and for models with a larger initial amount of angular momentum 
(MGEN2 and ACC2; Figure A2). 
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MGENOC; t=0 ACCO; t=0 



Figure Al. Same as Figure 3, but comparing the initial conditions of models MGENQc and ACC®. 
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MGEN2; t=0 ACC2; t=0 



Figure A2. Same as Figure 3, but comparing the initial conditions of models MGEN2 and ACC2. 
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APPENDIX B: SUPPLEMENTARY FIGURES 
(PHASE-SPACE MIXING) 

The additional figures presented in this section show, for the mod¬ 
els not illustrated in §4.1, the initial and final energy/z-angular mo¬ 
mentum distributions of the pristine and polluted populations. 
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Figure Bl. Same as Figures 6 and 7 but for model ACCS. 
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Figure B2. Same as Figures 6 and 7 but for model ACC2. 


Figure B3. Same as Figures 6 and 7 but for model MGENSa. 
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Figure B4. Same as Figures 6 and 7 but for model MGENSb. 


Figure B5. Same as Figures 6 and 7 but for model MGENQc. 
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Figure B6. Same as Figures 6 and 7 but for model MGEN2. 
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APPENDIX C: SUPPLEMENTARY FIGURES (VELOCITY 
DISPERSION) 

This section presents four additional figures to complement §4.2, 
Figure 10 and Figure 11. These illustrate the evolution of the x 
and z velocity dispersion profiles for models without net angular 
momentum initially (MGEN0C and ACCO; Figure Cl and Figure C3), 
and for models with a larger initial amount of angular momentum 
(MGEN2 and ACC2; Figure C2 and Figure C4). 
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Figure Cl. Same as Figure 10 but for models MGENSc and ACCO. 
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Figure C2. Same as Figure 10 but for models MGEN2 and ACC2. 
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Figure C3. Same as Figure 11 but for models MGENSc and ACC®. 
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Figure C4. Same as Figure 11 but for models MGEN2 and ACC2. 
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APPENDIX D: SUPPLEMENTARY FIGURES (ROTATION) 

This section presents two additional figures to complement §4.4 
and Figure 14 and illustrate the evolution of the rotation curve for 
models without net angular momentum initially (MGENOc and ACC0; 
Figure Dl), and for models with a larger initial amount of angular 
momentum (MGEN2 and ACC2; Figure D2). 
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Figure Dl. Same as Figure 14 but for models MGENSc and ACCO. 
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Figure D2. Same as Figure 14 but for models MGEN2 and ACC2. 
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